Delocalization and spreading in a nonlinear Stark ladder 
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We study the evolution of a wave packet in a nonlinear Schrodinger lattice equation subject 
to a dc bias. In the absence of nonlinearity all normal modes are spatially localized giving rise 
to a Stark ladder with an equidistant eigenvalue spectrum and Bloch oscillations. Nonlinearity 
induces frequency shifts and mode-mode interactions and destroys localization. With increasing 
strength of nonUnearity we observe: (I) localization as a transient, with subsequent subdiffusion 
(weak mode- mode interactions); (II) immediate subdiffusion (strong mode- mode interactions); (III) 
single site trapping as a transient, with subsequent explosive spreading, followed by subdiffusion. 
For single mode excitations and weak nonlinearities stability intervals are predicted and observed 
upon variation of the dc bias strength, which affect the short and long time dynamics. 
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The textbook solid state problem of an electron in a 
periodic potential with an additional electric field (see 
e.g. lead to investigations of Bloch oscillations Q 
and Landau- Zener tunneling [3] in differentphysical sys- 
tems. Optical waves in photonic lattices [4, 5] and ul- 
tracold atoms in optical lattices [6, 7] are recent exam- 
ples. If the well depth of the periodic potential is large 
enough, Landau-Zener tunneling is suppressed, the prob- 
lem is discretized (Stark ladder) and the resulting eigen- 
value problem is explicitly solved in terms of the local- 
ized eigenmodes [8*] of the system. Additional nonlinear 
terms in the wave equations are generated e.g. by a non- 
linear response of the medium (light) or by a mean-field 
treatment of quantum many-body interactions (BEG). 
These terms destroy the integrability of the wave equa- 
tions leading to decoherence and deterioration of Bloch 
oscillations, which was observed on time scales of a few 
oscillation periods [5, 7]. In the present work we study 
and explain the delocalization process induced by nonlin- 
ear terms. It leads to an incoherent subdiffusive spread- 
ing of any initially localized wave packet. 

We consider a discrete nonlinear Schrodinger equation 
(DNLS) with a dc bias E: 



(^n+l + ^n-l) + nE^n + /^l^nl'^n, (1) 



where the dot denotes the time derivative, is a com- 
plex amplitude of a wave at lattice site n, and P is the 
nonlinearity strength. Varying /3 is strictly equivalent to 
varying the norm Af = l^nP- We choose JV = 1. 

For /3 = Eq. ([1]) is reduced to an eigenvalue problem 
using the ansatz ^n(^) = exp(— At)^^. In the case of an 
infinite lattice the eigenvalues = Eu (with u being an 
integer) form an equidistant spectrum which extends over 
the whole real axis - the Stark ladder. The corresponding 
eigenfunctions or normal modes (NMs) obey the gener- 
alized translational invariance A, 



A^n^ i. With 

J^+i(x) + Jn-i(x) = {2n/x)Jn{x) for the Bessel function 



Jn{x) of the first kind [9] it follows Q 

4(0)^/ {-lTJ-n(2/E) , n<0 
^ Jn{2/E) , n > . 



(2) 



All NMs are spatially localized with an asymptotic de- 
cay lAnicx) I (l/£^) /n!, giving rise to the well-known 
Bloch oscillations with period Tb = 2tt/ E. In Fig. [1^ 
we plot the localization volume C = 1/X]^|^n^^|^ of 
an eigenstate, which characterizes the spatial extent of 
the NM, as a function of E . In the wide region of 
E e [0.003, 1] we find C = a-E'^ , a ^ 3.847 , b ^ 0.874. 
Note that for E = 10^ C is already close to its asymp- 
totic value £ = 1 for ^ ^ 00. The eigenvalue spacing E 
and the eigenvalue variation over a localization volume 
A = EjC determine the packet evolution in the presence 
of nonlinearity. 
With 



^n(t)=^(/>.(t)4^) 



Eq. ([T]) transforms into NM space 
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(4) 



(5) 



are the overlap integrals between the eigenfunctions ([2]). 
To characterize the wave packet dynamics we compute: 
i) the participation number P = 1/ l^'ul^ which mea- 
sures the number of strongest excited NMs; ii) the second 
moment m2 - Z^zy^l^z^P) which 

quantifies the degree of spreading of the wave packet and 
iii) the compactness index ( = P^/m2 which measures 
the degree of sparcity of the wavepacket. We use the 
SBAB2 symplectic integrator [10] to integrate Eq. ([T]) 
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FIG. 1: (Color online) (a) Black solid line: localization vol- 
ume C of the eigenfunction J^^^ versus E. Red dotted line: 
The fit of £ by a • with a = 3.847 and h = 0.874. In- 
sets: explicit form of the eigenfunction across the chain for 
two values ^ = 2 and ^ = 0.2 [(b), blue; (r), red]; (b) The 
diagram of the three regimes of spreading in the parameter 
space {f3,E). Empty and filled circles: numerically obtained 
thresholds which separate the three different regimes I, II, III. 
Lines connecting symbols guide the eye. Black dashed line: 
threshold between II and III obtained from the dimer model 
[see text]. In the limit of ^ ^ oo all lines merge to the 
asymptotic limit P oc E. 

numerically (see also [llj for details). The number of 
sites is chosen such as to exclude finite size effects in the 
evolution of the wave packets (typically from N = 500 to 
N = 1000). 

First, we study a single site initial excitation ^n(O) = 
6nO' In that case the amplitudes in NM space are ^zy(O) = 
Jy(2/E). The nonlinear frequency shift at site n = is 
S = p. We expect three qualitatively different regimes 
of spreading: (I) 5 < E', (II) E < 5 < A', (III) A < 
S. In case (I) the nonlinear frequency shift is less than 
the spacing between excited modes. Therefore no initial 
resonance overlap is expected, and the dynamics may 
evolve as the one for /3 = at least for long times. In case 
(II) resonance overlap happens, and the packet should 
evolve differently. For case (III) S tunes the excited site 
out of resonance with the neighbouring NMs. Resonances 
with more distant NMs are possible, but the overlap with 
these NMs is the weaker the further away they are ([2j). 
Therefore for long times the excited site may evolve as 
an independent oscillator (trapping). The diagram of the 
three regimes of spreading in the parameter space (/?, E) 
is shown in Fig. [TJd. 

Let us start with case (III). For E = 2 and (5 > 8.9, the 
single site excitation stays trapped up to times t = 10^ 
without significantly spreading into any other site of the 
lattice (violet curve in Fig|2]). Slightly lowering P we ob- 
serve that the excitation is trapped up to a some time 
Te which sensitively depends on /3 and changes by many 
orders of magnitude e.g. between 10^ to 10^ in the nar- 
row interval p e (8.05,8.9) for E = 2 (FigE]). For times 
t > Te dLii explosive and spatially asymmetric spreading 
is observed on a time scale of one Bloch period Tg. The 
packet spreads in the direction of NMs with larger eigen- 



FIG. 2: (color online) Single site excitation for E — 2. Second 
moment m2 and participation number P versus time in log- 
log plots for different values of p inside the interval where 
an explosive delocalisation of the trapped regime occurs: p = 
8.15, 8.25, 8.5 [(o) orange; (g) green; (r) red]. p = 8 [(b), blue]: 
intermediate regime. P = 8.9 [(v), violet]: trapped regime. 
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FIG. 3: |^n(t)|^ as a function of time (in logarithmic scale) 
and lattice site for ^ = 2, /3 = 8.15 (upper graph) and E = 
0.2, p — 5.2 (lower graph). The zooms into the explosion 
region are plotted on a linear time scale. Bloch oscillations 
are observed right after the explosion time Te- 

values, which provide the possibility of resonant energy 
transfer from the single site excitation due to its positive 
nonlinear frequency shift S. For about ten Bloch peri- 
ods Tb the packet shows Bloch oscillations, which then 
quickly decohere. Finally the packet spreads incoherently 
and subdiffusively (Figs j2|3]) . The explosion time Te is 
not monotonously changing with p (cf. Fig|2]), which in- 
dicates intermittency, i.e. the single site excitation can 
be closer or farer away from some regular structures in 
phase space. That distance may in turn control the value 
of Te' For E = 2 and P = 8 the packet spreads from 
scratch (blue curve in Fig|2]). 

In regime II the nonlinearity is not strong enough to 
trap a single site excitation, but the excited NMs are res- 
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FIG. 4: (color online) Single site excitation for E — 2. vfi2 and 
P versus time (log- log plots) for i) /3 = 2.1 [(b), blue]: weak 
nonlinear regime; ii) = 4 [(g), green]: intermediate regime 
and iii) /3 = 9 [(r), red]: trapped regime. Orange dashed line: 
7712 - Inset: C for /3 = 4. 



onantly interacting. Therefore the packet expands subd- 
iffusively from the very beginning (green curve in Fig. 2]). 
Since the eigenmodes of the linear system are spatially 
localized, the incoherent way of subdiffusive spreading re- 
sembles the one observed in disordered nonlinear systems 
lll,[l2|. The second moment 7772 appears to grow without 
bounds according to 7772 ~ . We fitted the exponent 
a for 20 different runs for £^ = 2 and = 4. Each run 
was obtained as follows. We integrated the equations up 
to t = 10^. The obtained distribution was transformed 
to NM space. The phases of all NMs were randomized, 
but the amplitudes were kept constant. Finally, we per- 
formed an inverse transformation to real space and con- 
tinued the integration up to 10^. As a result we obtained 
the average exponent a = 0.38 ±0.015. Let us reconsider 
(jl]) and keep only resonant terms v ^ v\ — V2 — = 
The resulting resonant normal form equations read 



(6) 



If at least two neighbouring NMs are excited, the reso- 
nant normal form (|6]) will have a connectivity similar to 
the original lattice equations, which allows spreading over 
the whole lattice. Excluding the possibility of integrabil- 
ity of dni, we expect this spreading to happen indeed, as 
observed in the numerical runs. Following [11] we con- 
clude, that all NMs are in resonance with their neigbour- 
hood. If the wave packet is in a state of equipartition 
at any time, it follows a = 0.5 [11]. In such an assumed 
equipartitioned state the compactness index ( ^ 3 [11]. 
However, the compactness index ( decreases in time way 
beyond the equipartition threshold (inset in Fig. [4j) and 
therefore the assumption of equipartitioned wave pack- 
ets is incorrect. This is the reason for the discrepancy 
between the expected (0.5) and observed (0.38) values of 
the exponent a. 

In regime I the nonlinear frequency shift 6 < E. No 
strong initial resonance overlap between NMs is expected. 
In that case the dynamics may (at least for long times) 



evolve as in the linear case. Eventually, the subdiffusive 
spreading takes place after some waiting time r (Fig. [4]). 
Note that the crossover from localization to spreading 
around time r is not accompanied by a rapid redistribu- 
tion of the energy between many neighboring sites as in 
the case of trapping with large nonlinear it ies. Therefore 
the onset of subdiffusion is simply related to the fact, that 
weakly chaotic dynamics is characterized by a small mo- 
mentary diffusion rate, and the corresponding diffusion 
process is observed only after times larger than r. 

In Fig. [TJd we plot the boundaries between the three 
regimes in the parameter space (/3, E). These boundaries 
are determined as follows. Keeping the value of E fixed 
and varying we integrate the equations up to t = 10^ 
starting from the single site excitation. Then we calculate 
7772 and P of the final wave packet. Being in regime I, 
these values are quite close to those for the linear case. 
For a critical value of /3 they start to grow indicating that 
we enter regime II. Finally, at a second and larger value 
of P the values of 1712 and P drop down again, signalling 
the entering of regime III. 

The border between regime II and III can be accounted 
for by using a dimer model. Since the trapped state cor- 
responds to only a single site being excited, we consider 
a dimer model where only one further lattice site (to the 
right) is added, in order to model the asymmetric energy 
transfer during the explosion. The resulting equations 
read 

which are solved using elliptic functions [13]. The dimer 
model is integrable, and the analytical computation of 
the initial condition |^n(0)P = ^n,o crossing a separatrix 
in the phase space of the dimer yields a relation between 
P and E. The transition to the trapped state at some 
critical nonlinearity P{E) is accompanied by divergence 
of the period of oscillations and above the threshold the 
energy is mostly concentrated on one site. The final equa- 
tion to be solved numerically reads 

-/?^£;+/3^(3£;^-i)+/3£;(2o-3£;^)+(£;^+4)^ = o. (7) 

As is seen from Fig(Tl3 it fits very well with the boundary 
between regime II and III obtained from the evaluation 
of the full set of equations. In particular, (|9]) correctly 
predicts that p{E ^ 0) ^ 4 and p{E oo) ^ E. 

A single mode excitation (j)iy{t = 0) = S^^o also exhibits 
the three different regimes of spreading. However, for 
small values of nonlinearity P a new intriguing feature 
of the short time dynamics follows. Indeed, consider the 
resonant normal form (|6]). For a single mode excitation, 
it follows strictly that no other NM is going to be excited. 
Therefore it is an exact solution of (|6]) 



0o(^) = exp (-i/5/o,o,o,o^) . 



(8) 



However, the full set of equations (j3|) will excite other 
NMs as well. These small perturbations may stay small 



or start to grow, depending on the stability of the solu- 
tion dH]) of equations ([6]). The linear equations for small 
perturbations with u read 

i^^y = P [/-i.,z.,o,o^o^-i. + 2/^,z/,o,o^i/] , (9) 

Note that only pairs of NMs with are coupled. With 

4>±j, ^ exp (-i/3/o,o,o,o^ + Jut) one gets 

= \P\ 7^-.,.,o,o-(2Wo-^o,o,o,o)'. (10) 

Thus if /?.^,i,,o,o > (2/zy,iy,o,o - ^o,o,o,o)^ it follows that 7^ 
is real, corresponding perturbations will grow exponen- 
tially fast in time, and the original single mode excitation 
is unstable. Since the overlap integrals depend on the dc 
bias E", we expect transitions from stability to instability 
upon varying E. Indeed, a numerical evaluation of that 
condition yields that the single mode excitation is stable 
in the limits: 1.134 < E < 1.23, 1.314 < E < 1.917, 
E > 2.314; otherwise it is unstable. In the instability 
window 1.917 < E < 2.314 (1.23 < E < 1.314) the 
modes u = ^1 {u = ±2) become unstable whereas the 
other modes remain stable. This behavior is confirmed 
by numerical simulations presented in Fig. [5] on short 
time scales. After sufficient exponential growth the am- 
plitude of the unstable mode saturates due to nonlinear 
corrections. Due to recurrences the energy eventually 
comes back to the initially excited mode after possible 
long times. We also find a correlation between the short 
time stability properties of a single mode excitation and 
the long time evolution of the wave packet. If the single 
mode excitation is launched in a stability window, the 
wave packet is practically not spreading up to t = 10^ 
(inset in Fig|5|). However, a small change of the dc bias 
value E can tune the system into an instability window. 
As a result, the wave packet shows subdiffusive spread- 
ing on the accessible time scales, starting with t ~ 10^ 
(inset in Figl5|). In conclusion, we have investigated var- 
ious regimes of wave packet spreading for the nonlinear 
Stark ladder. For large nonlinearities there is a long lived 
trapped regime with an explosive transition to Bloch os- 
cillations, followed by a subdiffusive spreading at large 
time scales. For intermediate nonlinearities we find subd- 
iffusion from scratch. For small nonlinearities the system 
is showing localization as a transient, with a subsequent 
subdiffusive spreading. For single mode excitations ad- 
ditional stability and instability intervals with respect to 
the dc bias strength exist. For weak nonlinearity this 
leads to an additional sensitive dependence of the tran- 
sient localization times on the value of the dc bias. In all 
cases, we observe that nonlinearity destroys integr ability, 
introduces chaos, and ultimately leads to a destruction 
of localization, which was a hallmark of the linear sys- 
tem due to phase coherence and destructive/constructive 
wave interference. 
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FIG. 5: (color online) Amplitudes of the excited, = 0, and 
neighboring modes, zy = 1 and z/ = 2, versus time for P = 0.4. 
(a) E=1.6, = [(b), blue], E=2, z/ = 0,1,2 [(o), orange; 
(r), red; (g), green], (b) E=1.18, zy - [(b), blue]; E=L26, 
zy = 0, 1,2 [(o), orange; (r), red; (g), green]. Inset: long time 
evolution. m2 versus time in log- log plots: E= 1.18, 1.26 [(b), 
blue; (o), orange]. 
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